Nonparametric graphical model for counts

Although multivariate count data are routinely collected in many application areas, there is surprisingly little work developing flexible models for characterizing their dependence structure. This is particularly true when interest focuses on inferring the conditional independence graph. In this article, we propose a new class of pairwise Markov random field-type models for the joint distribution of a multivariate count vector. By employing a novel type of transformation, we avoid restricting to non-negative dependence structures or inducing other restrictions through truncations. Taking a Bayesian approach to inference, we choose a Dirichlet process prior for the distribution of a random effect to induce great flexibility in the specification. An efficient Markov chain Monte Carlo (MCMC) algorithm is developed for posterior computation. We prove various theoretical properties, including posterior consistency, and show that our COunt Nonparametric Graphical Analysis (CONGA) approach has good performance relative to competitors in simulation studies. The methods are motivated by an application to neuron spike count data in mice.


Introduction
Graphical models provide an appealing framework to characterize dependence in multivariate data X i = (X i1 , . . ., X iP ) in an intuitive way.This article focuses on undirected graphical models or Markov random fields (MRFs).In this approach, each random variable is assigned as a node of a graph G which is characterized by the pair (V, E).Here V and E denote the set of nodes and set of connected edges of the graph G, with V = {1, . . ., P } and E ⊆ V × V .The graph G encodes conditional independence relationships in the data.We say X l and X k are conditionally independent if P (X l , X k |X −(l,k) ) = P (X l |X −(l,k) )P (X k |X −(l,k) ), with X −(l,k) denoting all random variables excluding X l and X k .Conditional independence between two random variables is equivalent to the absence of an edge between those two corresponding nodes in the graph.Thus the conditional independence of X l and X k would imply that the edge (k, l) is not present i.e. (k, l) / ∈ E.
Although there is a rich literature on graphical models, most of the focus has been specifically on Gaussian graphical models.For bounded discrete data, Ising (Ravikumar et al., 2010;Kolar et al., 2010) and multinomial graphical models (Jalali et al., 2011) have been studied.However, for unbounded count-valued data, the existing literature is limited.
Multivariate count data are routinely collected in genomics, sports, imaging analysis and text mining among many other areas, but most of the focus has been on latent factor and covariance structure models (Wedel et al., 2003;Zhou et al., 2012).The goal of this article is to address this gap and provide a flexible framework for statistical inference in count graphical models.
Besag first introduced pair-wise graphical models, deemed 'auto-models' in his seminal paper on MRFs (Besag, 1974).To define a joint distribution on a spatial lattice, he started with an exponential family representation of the marginal distributions and then added first order interaction terms.In the special case of count data, he introduced the Poisson auto-model.In this approach, the random variable at the i-th location X i follows a conditional Poisson distribution with mean µ i , dependent on the neighboring sites.Then µ i is given the form µ i = exp(α i + j β ij X j ).It can be shown that this conditional density model admits a joint density if and only if β ij ≤ 0 for all pairs of (i, j).Hence, only non-negative dependence can be accommodated.Gamma and exponential auto-models also have the same restriction due to non-negativity of the random variables.Yang et al. (2013) proposed to truncate the count support within the Poisson auto-model to allow both positive and negative dependence, effectively treating the data as ordered categorical.Allen and Liu (2012) proposed to fit the Poisson graphical model locally in a manner that allows both positive and negative dependencies, but this approach does not address the problem of global inference on G.
In the literature on spatial data analysis, many count-valued spatial processes have been proposed, but much of the focus has been on including spatial random effects instead of an explicit graphical structure.De Oliveira ( 2013) considered a random field on the mean function of a Poisson model to incorporate spatial dependence.However, conditional independence or dependence structure in the mean does not necessarily represent that of the data.The Poisson-Log normal distribution, introduced by Aitchison and Ho (1989), is popular for analyzing spatial count data (Chan and Ledolter, 1995;Diggle et al., 1998;Chib and Winkelmann, 2001;Hay and Pettitt, 2001).Here also the graph structure of the mean does not necessarily represent that of the given data.Hence, these models cannot be regarded as graphical models for count data.To study areal data, conditional autoregressive models (CAR) have been proposed (Gelfand and Vounatsou, 2003;De Oliveira, 2012;Wang and Kockelman, 2013).Although these models have an MRF-type structure, they assume the graph G is known based on the spatial adjacency structure, while our focus is on inferring unknown G.
Gaussian copula models are popular to characterize multivariate non-normal data (Xue-Kun Song, 2000;Murray et al., 2013).Mohammadi et al. (2017) developed a computational algorithm to build graphical models based on Gaussian copulas using methods developed by Dobra et al. (2011).However, it is difficult to model multivariate counts with zero inflated or multimodal marginals using a Gaussian copula.
Within a semiparametric framework, Liu et al. (2009) proposed a nonparanormal graphical model in which an unknown monotone function of the observed data follows a multivariate normal model with unknown mean and precision matrix subject to identifiability restrictions.This model has been popular for continuous data, providing a type of Gaussian copula.For discrete data the model is not directly appropriate, as mapping discrete to continuous data is problematic.To the best of our knowledge, there has been no work on nonparanormal graphical models for counts.In general conditional independence cannot be ensured if the function of the random variable is not continuous.For example if f is not monotone continuous, then conditional independence of X and Y does not ensure conditional independence of f (X) and f (Y ).
Apart from proposing a flexible graphical model for counts, we also aim at developing an efficient Bayesian computation scheme.Bayesian computation for Gaussian graphical models (GGMs) is somewhat well-developed (Dobra et al., 2011;Wang, 2012Wang, , 2015;;Mohammadi et al., 2015).Unfortunately, outside of GGMs likelihood-based inference is often problematic due to intractable normalizing constants in the likelihood functions.For example, the normalizing constant in the Ising model is too expensive to compute except for very small P .There are approaches related to surrogate likelihood (Kolar and Xing, 2008) by the relaxation of the log-partition function (Banerjee et al., 2008).Kolar et al. (2010) use conditional likelihood.Besag (1975) (Comets, 1992;Jensen and Künsch, 1994;Mase, 2000;Baddeley and Turner, 2000).Pensar et al. (2017) showed consistency of marginal pseudo-likelihood for discrete valued MRFs in a Bayesian framework.
Recently Dobra et al. ( 2018) used pseudo-likelihood for estimation of their Gaussian copula graphical model.Although pseudo-likelihood is popular in the frequentist domain for count data (Inouye et al., 2014;Ravikumar et al., 2010;Yang et al., 2013), its usage is still not very common in Bayesian estimation for count-valued MRFs.This is mainly because calculating conditional densities is expensive for count data due to unbounded support, making posterior computations hard to conduct.We implement an efficient Markov Chain Monte Carlo (MCMC) sampler for our model using pseudo-likelihood and pseudo-posterior formulations.Our approach relies on a provably accurate approximation to the normalizing constant in the conditional likelihood.We also provide a bound for the approximation error due to the evaluation of the normalizing constant numerically.
In Section 2, we introduce our novel graphical model.In Section 3, some desirable theoretical results are presented.Then we discuss computational strategies in Section 4 and present simulation results in Section 5. We apply our method to neuron spike data in mice in Section 6.We end with some concluding remarks in Section 7.

Modeling
Before introducing the model, we define some of the Markov properties related to the conditional independence of an undirected graph.A clique of a graph is the set of nodes where every two distinct nodes are adjacent i.e. connected by an edge.Let us define The joint distribution is locally Markov if for any j, X j and X l such that l ∈ V \ {N (j), j} are conditionally independent.If for three disjoint sets A, B and C of V , X A and X B are independent given X C whenever A and B are separated by C, the distribution is called globally Markov.The joint density is pair-wise Markov if for any i, j ∈ V such that (i, j) / ∈ E, X i and X j are conditionally independent.
We consider here a pair-wise MRF (Wainwright et al., 2007;Chen et al., 2014) which implies the following joint probability mass function (pmf) for the P dimensional random variable X, where f (X i ) is called a node potential function, f (X j , X l ) an edge potential function and we have f (X j , X l ) = 0 if there is no edge (j, l).Thus this distribution is pair-wise Markov Completing a specification of the MRF in (2.1) requires an explicit choice of the potential functions f (X j ) and f (X j , X l ).In the Gaussian case, one lets f (X j ) = −α j X 2 j and f (X j , X l ) = −β jl X j X l , where α j and β jl correspond to the diagonal and off-diagonal elements of the precision matrix Σ −1 = cov(X) −1 .In general, the node potential functions can be chosen to target specific univariate marginal densities.If the marginal distribution is Poisson, the appropriate node potential function is f (X j ) = α j X j − log(X j !).One can then choose the edge potential functions to avoid overly restrictive constraints on the dependence structure, such as only allowing non-negative correlations.Yang et al. (2013) identify edge potential functions with these properties for count data by truncating the support; for example, to the range observed in the sample.This reduces ability to generalize results, and in practice estimates are sensitive to the truncation level.We propose an alternative construction of the edge potentials that avoids truncation.

Model
We propose the following modified pmf for P -dimensional count-valued data X, ) using the notation of (2.1).
Lemma 1.Let F (•) be uniformly bounded by U , then the normalizing constant, say A(α, β), can be bounded as, By elementary calculations, these bounds can be obtained.The constant A(α, β) is the sum of the above pmf over the support of X.The sum reduces to a product of P many exponential series sums after replacing the function F (•) by its maximum.
Thus by modifying the edge potential function in this way using a bounded function of X, we can allow unrestricted support for all the parameters, allowing one to estimate both positive and negative dependence.Under the monotonicity restrictions on F (•), inference on the conditional independence structure tends to be robust to the specific form chosen.
We let F (•) = (tan −1 (•)) θ for some positive θ ∈ R + to define a flexible class of monotone increasing bounded functions.The parameter θ imposes additional flexibility and is calibrated to minimize differences in covariance between F (X) and X. Figure 1 illustrates how θ controls the range and shape of F (•). Figure 2 shows how the difference between covariances of F (X) and X vary for different values of θ in sparse and non sparse data cases.In both cases, the difference function has a unique minimizer.
Letting X t denote the t th independent realization of X, for t = 1, . . ., n, the pmf is where α tj 's are coefficients of different node potential functions and β jl 's are coefficients of the edge potential functions as before.We vary α tj with t to allow more flexibility in modeling the marginal densities.If β jl = 0, then X tj and X tl are conditionally independent for all t.We call our proposed method COunt Nonparametric Graphical Analysis (CONGA).Now we reparametrize (2.2) using log(λ tj ) = α tj and rewrite the model as, Pr(X t1 , . . ., X tP ) This reparametrizated model is more intuitive to understand.Due to the Poisson type marginal in (2.3), this model is suitable for data with over-dispersed marginals with respect to the Poisson at each node.Over-dispersion is typical in broad applications.We consider this reparametrized model in the rest of the paper.

Prior structure
To proceed with Bayesian computation, we put priors on the parameters.We have two set of parameters in (2.3), β and λ.For the β jl parameters, we choose simple iid Gaussian priors.It is straightforward to consider more elaborate shrinkage or variable selection priors for the β jl 's, but we find usual Gaussian priors have good performance in small to moderate-dimensional applications The parameter λ tj 's represent random effects; these parameters are not individually identifiable and are given random effects distributions λ tj ∼ D j .The distribution D j controls over-dispersion and the shape of the marginal count distribution for the j th node.To allow these marginals to be flexibly determined by the data, we take a Bayesian nonparametric approach using Dirichlet process priors D j ∼DP(M j D 0 ), with D 0 a Gamma base measure and M j a precision parameter, having M j ∼Ga(c, d) for increased data adaptivity.

Theoretical properties
We explore some of the theoretical properties of our proposed CONGA method.
Theorem 2. If we have β jl = 0, then X tj and X tl are conditionally independent for all t under (2.3).
This result is easy to verify by simply calculating the conditional probabilities.The details of the proof are in the Appendix.
Theorem 3. Any X that follows the probability mass function in (2.3) cannot be under dispersed with respect to the Poisson distribution.
We study posterior consistency under a fixed P and increasing n regime, assuming the prior of Section 2.2 with prespecified θ.Let the parameter space for G j be G j and that for β be R q , where q = P (P − 1)/2.Thus the complete parameter space We consider the prior Γj on G j and χ on β.
Let κ 0 be the truth for κ.We make following assumptions. Assumptions 2. For some large C > 0, let Q = {β : β ∞ < C}, where • ∞ stands for the infinity norm.Then β 0 ∈ Q and β 0 is in the support of χ.
We show that the truth belongs to the Kullback-Leibler support of the prior.Thus the posterior probability of any neighbourhood around the true p.m.f converges to one in P (n) κ 0 -probability as n goes to ∞ as a consequence of Schwartz (1965).Here P (n) κ is the distribution of a sample of n observations with parameter κ.Hence, the posterior is weakly consistent.The posterior is said to be strongly consistent if the posterior probability of any neighbourhood around the true p.m.f convergences to one almost-surely.Support of the data is a countable space.The weak and strong topologies on countable spaces are equivalent by Scheffe's theorem.In particular, weak topology and total variation topology are equivalent for discrete data.Weak consistency implies strong consistency.Thus the posterior for κ is also strongly consistent at κ 0 .A detail proof is in the Appendix.
Instead of assuming bounded support on the true distribution of random effects, one can also assume it to have sub-Gaussian tails.The posterior consistency result still holds with minor modifications in the current proof.

Computation
As motivated in Section 2.1, we estimate θ to minimize the differences in the sample covariance of X and F (X).In particular, the criteria is to minimize cov(tan −1 (X) θ ) − cov(X) F .This is a simple one dimensional optimization problem, which is easily solved numerically.
To update the other parameters, we use an MCMC algorithm, building on the approach of Roy et al. (2018).We generate proposals for Metropolis-Hastings (MH) using a Gibbs sampler derived under an approximated model.To avoid calculation of the global normalizing constant in the complete likelihood, we consider a pseudo-likelihood corresponding to a product of conditional likelihoods at each node.This requires calculations of P local normalizing constants which is computationally tractable.
The conditional likelihood at the j-th node is, We truncate this sum at a sufficiently large value B for the purpose of evaluating the conditional likelihood.The error in this approximation can be bounded by where CP (x, l) is the cumulative distribution function of the Poisson distribution with mean l evaluated at x.The above bound can in tern be bounded by a similar expression with (tan −1 (X tl )) θ replaced by (π/2) θ .One can tune B based on the resulting bound on the approximation error.In our simulation setting, even B = 70 makes the above bound numerically zero.We use B = 100 as a default choice for all of our computations.
We update λ .jusing the MCMC sampling scheme described in the Chapter 5 of Ghosal and Van der Vaart (2017) for the Dirichlet process mixture prior of λ ij based on the above conditional likelihood.For clarity this algorithm is described below: (i) First calculate the probability vector Q j for each j, such Q j (k) = Pois(X ij , λ kj ) and Otherwise we sample a new value as described below.
(iv) M j is sampled from Gamma(c + U, d − log(δ)), where U is the number of unique elements in λ .j, δ is sampled from Beta(M j , T ), and M j ∼Ga(c, d) a priori.
When we have to generate a new value for λ tj in step (iii), we consider the following scheme.
(i) Generate a candidate λ c tj from Gamma(a + X tj , b + 1).
(ii) Adjust the update λ c tj = λ 0 tj + K 1 (λ c tj − λ 0 tj ), where λ 0 tj is the current value and K 1 < 1 is a tuning parameter, adjusted with respect to the acceptance rate of the resulting Metropolis-Hastings (MH) step.
(iii) We use the pseudo-likelihood based on the conditional likelihoods in (4.1) to calculate the MH acceptance probability.
Thus s is the P × P gram matrix of (tan −1 X) θ , standardized over columns.
(i) Generate an update for Ω l,−l using the posterior distribution as in Wang (2012).Thus a candidate Ω c l,−l is generated from MVN(−Cs l,−l , C), where C = (( , where Ω 0 l,−l is the current value and K 2 is a tuning parameter, adjusted with respect to the acceptance rate of the following MH step.Also K 2 should always be less than (Ω c l,−l − Ω 0 l,−l ) 2 .
(iii) Use the pseudo-likelihood based on the conditional likelihoods in (4.1), multiplying over t to calculate MH acceptance probability.

Simulation
We consider four different techniques for generating multivariate count data.One approach is based on a Gaussian copula type setting.The other three are based on competing methods.We compare the methods based on false positive and false negative proportions.
We include an edge in the graph between the j th and l th nodes if the 95% credible interval for β jl does not include zero.There is a decision theoretic proof to justify such an approach in Thulin (2014).We compare our method CONGA with TPGM, SPGM, LPGM, huge, BDgraph and ssgraph.The first three are available in R package XMRF and the last two are in R packages BDgraph and ssgraph respectively.The function huge is from R package huge which fits a nonparanormal graphical model.The last two methods fit graphical models using Gaussian copulas and ssgraph uses spike and slab priors in estimating the edges.
To simulate data under the first scheme, we follow the steps given below.
(i) Generate n many multivariate normals of length c from MVN(0 c , Ω −1 c×c ), where 0 c is the vector of zeros of length c.This produces a matrix X of dimension n × c.
(ii) We calculate the matrix P n×c , which is P ij = Φ(X ij ), where Φ is the cumulative density function of the standard normal.
(iii) The Poisson random variable Y n×c is Y ij = QP (P ij , λ) for a given mean parameter λ with QP the quantile function of Poisson(λ). on modeling assumptions that CONGA avoids.CONGA beats BDgraph and ssgraph in almost all the scenarios in terms of false positive proportions.The false negative proportions are comparable.The function 'huge' performs similarly to CONGA when the data are generated using TPGM, SPGM and LPGM.But CONGA is better than all other methods when the data are generated using the Gaussian copula type setting.This is likely due to the fact that the other cases correspond to simulating data from one of the competitors models.

Neuron spike count application
The dataset records neuron spike counts in mice across 37 neurons in the brain under the influence of three different external stimuli, 2-D sinusoids with vertical gradient, horizontal gradient, and the sum.These neurons are from the same depth of the visual cortex of a mouse.The data are collected for around 400 time points.In Figure 3, we plot the marginal densities of the spike counts of four neurons under the influence of stimuli 0. We see that there are many variations in the marginal densities, and the densities are multi-modal for some of the cases.Marginally at each node, we also have that the variance is more than the corresponding mean for each of the three stimuli.

Estimation
We apply exactly the same computational approach as used in the simulation studies.To additionally obtain a summary of the weight of evidence of an edge between nodes j and l, we calculate S jl = |0.5 − P (β jl > 0)| /0.5, with P (β jl > 0) the posterior probability estimated from the MCMC samples.We plot the estimated graph with edge thickness proportional to the values of S jl .Thus thicker edges suggest greater evidence of an edge in Figures 4 to 6.To test for similarity in the graph across stimuli, we estimate 95% credible regions for ∆ s,s jl = β s jl − β s jl , denoting the difference in the (j, l) edge weight parameter under stimuli s and s , respectively.We flag those edges (j, l) having 95% credible intervals for ∆ s,s jl not including zero as significantly different across stimuli.

Inference
We find that there are 129, 199 and 110 connected edges respectively for stimuli 0, 1 and 2.
Among these edges, 38 are common for stimulus 0 and 1.The number is 15 for stimulus 0 and 2, and 28 for stimulus 1 and 2. There are 6 edges that are common for all of the stimuli.
We also list the most significant 10 edges for each stimulus in Table 5.We find that the node 27 is present in all of them.This node seems to have a significant interconnections network across stimulus.Here we find 82.13%similarity between the estimated weighted networks under the influence of stimulus 0 and 1.It is 84.98% for the pair 0 and 2. For 1 and 2, it is 79.43%.Stimulus 0 is a combination of stimuli 1 and 2. This could be the reason that the estimated graph under influence of stimulus 0 has the highest similarity with the other estimated graphs.
Figure 4: Estimated weighted network under the influence of stimuli 0. The edge width is proportional to the degree of significance.

Discussion
Our count nonparametric graphical analysis (CONGA) method is useful for multivariate count data, and represents a starting point for more elaborate models and other research di-   rections.One important direction is to time series data.In this respect, a simple extension is to define an autoregressive process for the baseline parameters α tj , inducing correlation in α t−1,j and α tj , while leaving the graph as fixed in time.A more elaborate extension would instead allow the graph to evolve dynamically by replacing the β jl parameters with β tjl , again defining an appropriate autoregressive process.
An additional interesting direction is to flexible graphical modeling of continuous positivevalued multivariate data.Such a modification is conceptually straightforward by changing the term log(X tj !) to the corresponding term in the gamma distribution.

Proof of Theorem 4
For q, q * ∈ the space of probability measure P, let the Kullback-Leibler divergences be given by K(q * , q) = q * log q * q V (q * , q) = q * log 2 q * q .
Thus the posterior is weakly consistent.The weak and strong topologies on countable spaces are equivalent by Scheffe's theorem.Thus the posterior for κ is also strongly consistent at κ 0 .
by construction.Then (2.1) satisfies the Hammersley-Clifford theorem(Hammersley and   Clifford, 1971), which states that a probability distribution having a strictly positive density satisfies a Markov property with respect to the undirected graph G if and only if its density can be factorized over the cliques of the graph.Since our pair-wise MRF is pair-wise Markov, we can represent the joint probability mass function as a product of mass functions of the cliques of graph G.The existence of such a factorization implies that this distribution has both global and local Markov properties.

Figure 3 :
Figure 3: Marginal densities of spike count of the four selected neurons under the influence of stimuli 0.

Figure 5 :
Figure 5: Estimated weighted network under the influence of stimuli 1.The edge width is proportional to the degree of significance.

Figure 6 :
Figure 6: Estimated weighted network under the influence of stimuli 2. The edge width is proportional to the degree of significance.

Table 1 :
Performance of the competing methods against our proposed method with 10 nodes.Top row indicates the method used to estimate and the first column indicates the method used to generate the data.p 1 and p 2 stand for false positive and false negative proportions.

Table 2 :
Performance of the competing methods against our proposed method with 30 nodes.Top row indicates the method used to estimate and the first column indicates the method used to generate the data.p 1 and p 2 stand for false positive and false negative proportions.

Table 3 :
Performance of the competing methods against our proposed method with 50 nodes.Top row indicates the method used to estimate and the first column indicates the method used to generate the data.p 1 and p 2 stand for false positive and false negative proportions.

Table 4 :
Top 5 nodes with maximum number of connected edges under the influence of stimuli 0, 1 and 2 are listed below.

Table 5 :
Top 10 most significant edges under the influence of stimulus 0, 1 and 2 with 1 as the estimated measure of significance are listed below.